All runs Ulva PHOTOSYNTHESIS Analysis, Script Chunks, and Plots

This is the analysis of all the Ulva lactuca salinity and nutrient experiments conducted on the lanai in St. John 616 from September 2021 to October 2022. These experiments incorporated four paired salinity and nutrient treatments with three temperatures. Each of the first four runs produced an n = 2 and was repeated initially 8 times for a total of n = 16. Data gaps were identified and filled in February, April, and October 2022. This output reflects all data totaling six treatments for Ulva lactuca.

Packages loaded:

library(lme4)
library(lmerTest)
library(effects)
library(car)
library(MuMIn)
library (dplyr)
library(emmeans)
library(DHARMa)
library(performance)
library(patchwork)
library(rstatix)
#for plots and tables
library(ggplot2)
library(ggpubr)
library(forcats)
library(RColorBrewer)
library(tidyverse)
library(sjPlot)
library(sjmisc)
#library(mmtable2)
library(gt)
library(purrr)
library(stringr)
library(tidyr)

Load and prepare the dataset

Open the output dataset generated by the ps_script_clean_to_ek_alpha.R script in the phytotools_alpha_ek project This file was normalized to quantum efficiency of photosynthesis as this seems to be more accurate (changed in fitWebb input) per Silsbe and Kromkamp (2012)

all_runs_photosyn_data <- read.csv("../data_input/hyp_ulva_all_runs_ek_alpha_normalized.csv")

Assign run as a factor

all_runs_photosyn_data$Run <- as.factor(all_runs_photosyn_data$Run)

Assign temperature as a factor

all_runs_photosyn_data$Temperature <- as.factor(all_runs_photosyn_data$Temp...C.)

Assign treatment as characters from integers then to factors

all_runs_photosyn_data$Treatment <- as.factor(as.character(all_runs_photosyn_data$Treatment))

Assign deltaNPQ as a factor

all_runs_photosyn_data$deltaNPQ <- as.factor(all_runs_photosyn_data$deltaNPQ)

Subset the data and toggle between the species for output. Use Day 9 for final analysis ONLY This will also assign the proper labels for plots

ulva <- subset(all_runs_photosyn_data, Species == "ul" & RLC.Day == 9)
ulva$treatment_graph[ulva$Treatment == 0] <- "1) 35ppt/0.5umol"
ulva$treatment_graph[ulva$Treatment == 1] <- "2) 35ppt/14umol" 
ulva$treatment_graph[ulva$Treatment == 2] <- "3) 28ppt/27umol" 
ulva$treatment_graph[ulva$Treatment == 3] <- "5) 18ppt/53umol" 
ulva$treatment_graph[ulva$Treatment == 4] <- "6) 11ppt/80umol"
ulva$treatment_graph[ulva$Treatment == 2.5] <- "4) 28ppt/53umol"

Add a column for growth rate from growth rate dataset to the already subsetted Ulva data frame

growth_rate <- read.csv("/Users/Angela/src/work/limu/algal_growth_photosynthesis/data_input/all_runs_growth_102222.csv")
growth_rate$Species <- as.factor(growth_rate$Species)
growth_rate$treatment <- as.factor(growth_rate$treatment)

Subset for Ulva only and calculate growth rate from final and initial weights

gr_ulva <- subset(growth_rate, Species == "Ul")
ulva$growth_rate <- round((gr_ulva$final.weight - gr_ulva$Initial.weight) / gr_ulva$Initial.weight * 100, digits = 2)

#Run the model Run model for rETRmax with two fixed effect variables and three random effects variables

#run model without interaction between the treatments and temperature
all_runs_photosyn_model_noint <- lmer(formula = rETRmax ~ Treatment + Temperature + (1 | Run)
                                      + (1 | Plant.ID) + (1 | RLC.Order), data = ulva)

rETRmax – Make a histogram and residual plots of the data

hist(ulva$rETRmax, main = paste("Ulva lactuca rETRmax"), col = "olivedrab3", labels = TRUE)

#or
ulva %>% ggplot(aes(rETRmax)) +
  geom_histogram(binwidth=5, fill = "#5BB300", color = "black", size = 0.25, alpha = 0.85) +
  theme_bw()

plot(resid(all_runs_photosyn_model_noint) ~ fitted(all_runs_photosyn_model_noint))

qqnorm(resid(all_runs_photosyn_model_noint))
qqline(resid(all_runs_photosyn_model_noint))

rETRmax – Check the performance of the model

performance::check_model(all_runs_photosyn_model_noint)

These outputs show the model is acceptable

rETRmax – Check r2 for model fit and print the model statistics summary

r.squaredGLMM(all_runs_photosyn_model_noint)
## Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help page.
##            R2m       R2c
## [1,] 0.6072108 0.7014082
summary(all_runs_photosyn_model_noint)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: rETRmax ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) +  
##     (1 | RLC.Order)
##    Data: ulva
## 
## REML criterion at convergence: 2190.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8257 -0.6193  0.1037  0.5129  3.4274 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  Plant.ID  (Intercept)  20.433   4.520  
##  Run       (Intercept)   7.118   2.668  
##  RLC.Order (Intercept)   6.907   2.628  
##  Residual              109.224  10.451  
## Number of obs: 288, groups:  Plant.ID, 143; Run, 8; RLC.Order, 6
## 
## Fixed effects:
##               Estimate Std. Error     df t value Pr(>|t|)    
## (Intercept)     35.172      3.016  9.396  11.661 6.72e-07 ***
## Treatment1      20.950      3.234  7.238   6.478 0.000296 ***
## Treatment2      22.035      3.234  7.238   6.814 0.000215 ***
## Treatment2.5    42.294      4.010  3.766  10.547 0.000622 ***
## Treatment3      38.157      3.234  7.238  11.799 5.48e-06 ***
## Treatment4      39.481      3.234  7.238  12.208 4.33e-06 ***
## Temperature27   -4.863      2.400 29.820  -2.027 0.051735 .  
## Temperature30   -2.783      2.136 67.783  -1.303 0.197127    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Trtmn1 Trtmn2 Trt2.5 Trtmn3 Trtmn4 Tmpr27
## Treatment1  -0.642                                          
## Treatment2  -0.642  0.782                                   
## Treatmnt2.5 -0.518  0.483  0.483                            
## Treatment3  -0.642  0.782  0.782  0.483                     
## Treatment4  -0.642  0.782  0.782  0.483  0.782              
## Temperatr27 -0.377  0.002  0.002  0.000  0.002  0.002       
## Temperatr30 -0.362  0.000  0.000  0.000  0.000  0.000  0.473

Use Bartlett’s test to check for equal variance

bartlett.test(rETRmax ~ Treatment, data = ulva)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  rETRmax by Treatment
## Bartlett's K-squared = 12.584, df = 5, p-value = 0.02761

Run Welch’s ANOVA if not equal variances

welch_anova_treatment <- oneway.test(rETRmax ~ Treatment, data = ulva, var.equal = FALSE)
welch_anova_treatment
## 
##  One-way analysis of means (not assuming equal variances)
## 
## data:  rETRmax and Treatment
## F = 128.4, num df = 5.00, denom df = 130.81, p-value < 2.2e-16
welch_anova_temp <- oneway.test(rETRmax ~ Temperature, data = ulva, var.equal = FALSE)
welch_anova_temp
## 
##  One-way analysis of means (not assuming equal variances)
## 
## data:  rETRmax and Temperature
## F = 1.4155, num df = 2.00, denom df = 189.85, p-value = 0.2454
games_howell_test(ulva, rETRmax ~ Treatment, conf.level = 0.95, detailed = TRUE)

rETRmax – Plot and make a table of the results. Also get means for the treatments

plot(allEffects(all_runs_photosyn_model_noint))

tab_model(all_runs_photosyn_model_noint)
  rETRmax
Predictors Estimates CI p
(Intercept) 35.17 29.23 – 41.11 <0.001
Treatment [1] 20.95 14.58 – 27.32 <0.001
Treatment [2] 22.04 15.67 – 28.40 <0.001
Treatment [2.5] 42.29 34.40 – 50.19 <0.001
Treatment [3] 38.16 31.79 – 44.52 <0.001
Treatment [4] 39.48 33.11 – 45.85 <0.001
Temperature [27] -4.86 -9.59 – -0.14 0.044
Temperature [30] -2.78 -6.99 – 1.42 0.194
Random Effects
σ2 109.22
Ï„00 Plant.ID 20.43
Ï„00 Run 7.12
Ï„00 RLC.Order 6.91
ICC 0.24
N Run 8
N Plant.ID 143
N RLC.Order 6
Observations 288
Marginal R2 / Conditional R2 0.607 / 0.701 
ulva %>% group_by(Treatment) %>% summarise_at(vars(rETRmax), list(mean = mean))
ulva %>% ggplot(aes(treatment_graph, rETRmax)) + 
  geom_boxplot(size=0.5) + 
  geom_point(alpha = 0.5, size = 3, aes(color = Temperature), show.legend = FALSE) + 
  labs(x="salinity/nitrate", y= "rETRmax (μmols electrons m-2 s-1)", title= "A", subtitle = "Ulva lactuca") + 
  scale_x_discrete(labels = c("35ppt/0.5umolN", "35ppt/14umolN", "28ppt/27umolN", "28ppt/53umolN", "18ppt/53umolN", "11ppt/80umolN")) + 
  ylim(-1, 160) + stat_mean() + 
  geom_hline(yintercept=0, color = "purple", size = 0.5, alpha = 0.5) +
  theme_bw() +
  theme(plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05), plot.subtitle = element_text(face = "italic", vjust = -20, hjust = 0.05))

Plot a regression between the photosynthetic independent variables of interest and growth rate

ulva_growth_etr_graph <- ggplot(ulva, aes(x=rETRmax, y=growth_rate)) + 
  geom_point(alpha = 0.5, size = 3, show.legend = TRUE, aes(color = Treatment)) + 
  geom_smooth(method = "lm", col = "black") + theme_bw() + 
  labs(title = "Ulva lactuca rETRmax vs Growth Rate", x = "rETRmax (μmols electrons m-2 s-1)", 
       y = "growth rate (%)") + stat_regline_equation(label.x = 25, label.y = 165) + stat_cor()
ulva_growth_etr_graph
## `geom_smooth()` using formula 'y ~ x'

Temperature did not have a significant effect on the outcome for rETRmax

Ek – Run the model

Run model for minimum saturating irradiance (Ek) with two fixed effect variables and three random effects variables

all_runs_photosyn_model_ek <- lmer(formula = ek.1 ~ Treatment + Temperature + (1 | Run)
                                      + (1 | Plant.ID) + (1 | RLC.Order), data = ulva)

Ek – Make a histogram and residual plots of the data for ulva

hist(ulva$ek.1, main = paste("Ulva lactuca Ek"), col = "darkolivegreen3", labels = TRUE)

plot(resid(all_runs_photosyn_model_ek) ~ fitted(all_runs_photosyn_model_ek))

qqnorm(resid(all_runs_photosyn_model_ek))
qqline(resid(all_runs_photosyn_model_ek))

ulva %>% ggplot(aes(ek.1)) +
  geom_histogram(binwidth=5, fill = "#5BB300", color = "black", size = 0.25, alpha = 0.85) +
  theme_bw()

Ek – Check the performance of the model

performance::check_model(all_runs_photosyn_model_ek)

Ek – Check r2 for model fit and print the model statistics summary

r.squaredGLMM(all_runs_photosyn_model_ek)
##            R2m      R2c
## [1,] 0.5492483 0.657886
summary(all_runs_photosyn_model_ek)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: ek.1 ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) +  
##     (1 | RLC.Order)
##    Data: ulva
## 
## REML criterion at convergence: 2296.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.5526 -0.5975  0.1036  0.5529  4.3552 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  Plant.ID  (Intercept)  31.629   5.624  
##  Run       (Intercept)   5.569   2.360  
##  RLC.Order (Intercept)  13.195   3.632  
##  Residual              158.693  12.597  
## Number of obs: 288, groups:  Plant.ID, 143; Run, 8; RLC.Order, 6
## 
## Fixed effects:
##               Estimate Std. Error     df t value Pr(>|t|)    
## (Intercept)     50.512      3.400 11.858  14.856 5.04e-09 ***
## Treatment1      22.632      3.451  9.064   6.557 0.000101 ***
## Treatment2      20.823      3.451  9.064   6.033 0.000189 ***
## Treatment2.5    40.137      4.035  3.766   9.947 0.000771 ***
## Treatment3      41.896      3.451  9.064  12.139 6.54e-07 ***
## Treatment4      43.942      3.451  9.064  12.732 4.34e-07 ***
## Temperature27   -7.842      3.010 34.223  -2.606 0.013480 *  
## Temperature30   -5.489      2.641 71.438  -2.078 0.041278 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Trtmn1 Trtmn2 Trt2.5 Trtmn3 Trtmn4 Tmpr27
## Treatment1  -0.576                                          
## Treatment2  -0.576  0.722                                   
## Treatmnt2.5 -0.492  0.485  0.485                            
## Treatment3  -0.576  0.722  0.722  0.485                     
## Treatment4  -0.576  0.722  0.722  0.485  0.722              
## Temperatr27 -0.417  0.002  0.002  0.000  0.002  0.002       
## Temperatr30 -0.398  0.000  0.000  0.000  0.000  0.000  0.470

Ek – Run Bartlett’s test and Welch’s ANOVA with Games Howell test for pairwise comparisons

bartlett.test(ek.1 ~ Treatment, data = ulva)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  ek.1 by Treatment
## Bartlett's K-squared = 12.071, df = 5, p-value = 0.03383
welch_anova_treatment <- oneway.test(ek.1 ~ Treatment, data = ulva, var.equal = FALSE)
welch_anova_treatment
## 
##  One-way analysis of means (not assuming equal variances)
## 
## data:  ek.1 and Treatment
## F = 89.233, num df = 5.00, denom df = 130.86, p-value < 2.2e-16
welch_anova_temp <- oneway.test(ek.1 ~ Temperature, data = ulva, var.equal = FALSE)
welch_anova_temp
## 
##  One-way analysis of means (not assuming equal variances)
## 
## data:  ek.1 and Temperature
## F = 3.2615, num df = 2.00, denom df = 189.43, p-value = 0.04049
games_howell_test(ulva, ek.1 ~ Treatment, conf.level = 0.95, detailed = TRUE)

Ek – Plots and tables for the results

tab_model(all_runs_photosyn_model_ek)
  ek.1
Predictors Estimates CI p
(Intercept) 50.51 43.82 – 57.21 <0.001
Treatment [1] 22.63 15.84 – 29.43 <0.001
Treatment [2] 20.82 14.03 – 27.62 <0.001
Treatment [2.5] 40.14 32.19 – 48.08 <0.001
Treatment [3] 41.90 35.10 – 48.69 <0.001
Treatment [4] 43.94 37.15 – 50.74 <0.001
Temperature [27] -7.84 -13.77 – -1.92 0.010
Temperature [30] -5.49 -10.69 – -0.29 0.039
Random Effects
σ2 158.69
Ï„00 Plant.ID 31.63
Ï„00 Run 5.57
Ï„00 RLC.Order 13.19
ICC 0.24
N Run 8
N Plant.ID 143
N RLC.Order 6
Observations 288
Marginal R2 / Conditional R2 0.549 / 0.658 
plot(allEffects(all_runs_photosyn_model_ek))

ulva %>% group_by(Treatment) %>% summarise_at(vars(ek.1), list(mean = mean))
ulva %>% ggplot(aes(treatment_graph, ek.1)) + 
  geom_boxplot(size=0.5) + 
  geom_point(alpha = 0.5, size = 3, aes(color = Temperature), show.legend = FALSE) + 
  labs(x="salinity/nitrate", y= "Ek (μmols photons m-2 s-1)", title= "A", subtitle = "Ulva lactuca") + 
  scale_x_discrete(labels = c("35ppt/0.5umolN", "35ppt/14umolN", "28ppt/27umolN", "28ppt/53umolN", "18ppt/53umolN", "11ppt/80umolN")) + 
  ylim(-1, 160) + stat_mean() + 
  geom_hline(yintercept=0, color = "purple", size = 0.5, alpha = 0.5) +
  theme_bw() +
  theme(plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05), plot.subtitle = element_text(face = "italic", vjust = -20, hjust = 0.05))
## Warning: Removed 1 rows containing non-finite values (stat_boxplot).
## Warning: Removed 1 rows containing non-finite values (stat_mean).
## Warning: Removed 1 rows containing missing values (geom_point).

plot linear regression Ek with growth rate

ulva_growth_etr_graph <- ggplot(ulva, aes(x=ek.1, y=growth_rate)) + 
  geom_point(alpha = 0.5, size = 3, show.legend = TRUE, aes(color = Treatment)) + 
  geom_smooth(method = "lm", col = "black") + theme_bw() + 
  labs(title = "Ulva lactuca Ek vs Growth Rate", x = "Ek (μmols photons m-2 s-1)", 
       y = "growth rate (%)") + stat_regline_equation(label.x = 25, label.y = 165) + stat_cor()
ulva_growth_etr_graph
## `geom_smooth()` using formula 'y ~ x'

Plot a regression between the photosynthetic independent variables of interest and growth rate

ulva_growth_etr_ek_graph <- ggplot(ulva, aes(x=rETRmax, y=ek.1)) + 
  geom_point(alpha = 0.5, size = 3, show.legend = TRUE, aes(color = Treatment)) + 
  geom_smooth(method = "lm", col = "black") + theme_bw() + 
  labs(title = "Ulva lactuca rETRmax vs Ek", x = "rETRmax (μmols electrons m-2 s-1)", 
       y = "Ek (μmols photons m-2 s-1)") + stat_regline_equation(label.x = 25, label.y = 165) + stat_cor()
ulva_growth_etr_ek_graph
## `geom_smooth()` using formula 'y ~ x'